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A quantitative comparison between experimental and Monte Carlo simulation results for the epi- 
taxial growth of Cu/Cu(001) in the submonolayer regime is presented. The simulations take into ac- 
count a complete set of hopping processes whose activation energies are derived from semi-empirical 
calculations using the embedded-atom method. The island separation is measured as a function of 
the incoming flux and the temperature. A good quantitative agreement between the experiment and 
simulation is found for the island separation, the activation energies for the dominant processes, and 
the exponents that characterize the growth. The simulation results are then analyzed at lower cover- 
ages, which are not accessible experimentally, providing good agreement with theoretical predictions 
as well. 
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The growth of a newly deposited layer (adlayer) of 
a thin metal film in molecular-beam epitaxy (MBE) 
through island nucleation involves three kinetic regimes 
with respect to the coverage 6: (a) Nucleation regime, 
dominated by island nucleation; (b) Aggregation regime, 
where newly deposited atoms (adatoms) are mostly cap- 
tured by existing islands. In this regime the system is in 
a quasi steady-state, a property that is employed in the 
scaling theory; (c) Coalescence regime, where islands are 
large enough to interact and merge. 

Experiments using scanning tunneling -microscopy 
(STM)ETj low-energy electron diffractioncm and other 
methods,u made in attempt to study the aggregation 
regime, reveal an exponential dependencC|-|iplj-.tiic island 
density N, on the inverse temperature 1/H Ml as well 
as a power law dependence on the flux i^.BQtra This be- 
havior can be described by the following formula:!! 



FY 

— exp (E cS /k B T) . 



(1) 



where I is the mean island separation related to N 
through N ~ £~ 2 . In the experimental-data analysis, 
7, v and E c g are fitting parameters, where E c r is a mea- 
sure of the activation energy for the rate-limiting move, 
v is the attempt rate for that move and 7 is a scaling ex- 
ponent. The theoretical challenge is, therefore, to relate 
7, v and E c f[ to the fundamental processes at the atomic 
level. 

In this paper we show that kinetic Monte Carlo (MC) 
simulations based on energy barriers obtained via the 
embedded-atom method (EAM), provide a good quan- 
titative description of diffusion and growth on Cu(001). 



To this end we performed kinetic MC simulations of is- 
land growth under conditions identical to those employed 
experimentally in Refs. In these experiments spot- 
profile analysis of low-energy electron diffraction (SPA- 
LEED) was used to measure the mean separation £ be- 
tween Cu islands on Cu(001). The measurements were 
taken as a function of the flux F for three temperatures, 
T = 213, 223 and 263 K (Fig. g), and as a function of the 
temperature T at a constant flux F = 3.21 x 10~ 4 ML ■ 
s" 1 (Fig. |). 

In our kinetic MC simulations atoms are randomly de- 
posited on a square-lattice substrate of 250 x 250 sites 
that corresponds to a terrace width of ~ 64 nrnJla These 
atoms attach irreversibly to the surface and hop as ran- 
dom walkers to unoccupied nearest-neighbor (NN) sites. 
Each hop involves an activation energy-barrier, E n , that 
depends on the configuration of occupied and unoccu- 
pied adjacent .sites in the 3x3 square around the hop- 
ping atom.li2ni3 When an atom is deposited on top of 
an island it is incorporated at a random position along 
the island perimeter. The nucleation of a second layer is 
thus suppressed, which is a good approximation for small 
island sizes at low coverages. The hopping-rate (in units 
of hops per second), h n , for some configuration n isllil 



exp(-E n /K B T), 



(2) 



where v is the attempt rate common to all moves and E n 
is calculated using the EAmMBB 

We have reconstructed the experimental conditions in 
our kinetic MC simulations and obtained for each exper- 
imental curve in Figs, [l] and || a corresponding simulated 
curve. Each simulated data point is an average over 20 
runs. In Fig. [j] we present the experimental results and 
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T(K) 


experiment 


simulations 


simulations 




9 = 0.3 (ML) 


9 = 0.3 (ML) 


9 = 0.125 (ML) 


213 


0.28 ± 0.02 


0.30 ± 0.01 


0.309 ± 0.006 


223 


0.42 ± 0.02 


0.326 ± 0.006 


0.319 ±0.006 


263 


0.54 ±0.01 


0.51 ± 0.01 


0.396 ± 0.006 



TABLE I. The values for the exponent 7 defined in Eq. (|l| 
for 9 = 0.3 ML (Fig. 0) and 9 = 0.125 ML. 



the corresponding simulation results, for the island sepa- 
ration I vs. 1/F at T = 213, 223 and 263 K. All data 
points were taken at 9 = 0.3 MLO Clearly, there is 
a good agreement between the experiment and simula- 
tion. The curves follow a power-law behavior according 
to Eq. (Q) in a range of one and -a. half decades, from 
which the exponent 7 is extracted. E£ -The results for the 
exponent 7 are presented in Table It is found to in- 
crease in the range 7 = 0.3 - 0.5 as the temperature is 
raised, suggesting that the system undergoes a crossver in 
its kinetic behavior.E3 Based on scaling theory, one may 
conclude from the results 7 ~ 1/3 that dimers and larger 
islands are stable and immobile, and adatoms are the only 
mobile entities on the surfaced However, care should be 
exercised when drawing such conclusions. Specifically, in 
our case the value 7 ~ 1 /3 is found in spite of significant 
dinner mobility as we show later. 

Unlike the energy barriers, we do not have a value for 
the attempt rate v, from atomic scale calculations. In- 
stead, the experimental value of v is obtained by fitting 
the simulation and experimental results. This is possible 
since in the simulation v sets a fundamental clock rate, 
while the simulation results depend only on the ratio v/F 
rather then on F and v separately. Using this property 
we perform the simulations for a broad range of values 
of v/F and plot the island separation £ vs. v/F. The 
attempt rate v, is then obtained as the value for which 
the three simulated curves in Fig. |l|, overlap simultane- 
ously the three experimental curves. It is found to be 
i/ = 1.2x 10 13 s" 1 . 

In Fig. H we present the experimental results (■) for the 
island separation I vs. 1/T, at F = 3.21 x ICL1 4 ML-s -1 
and the corresponding simulation results (*),til taken at 
9 = 0.3 ML. The apparent good agreement between them 
is a further confirmation of the correctness of v. Both 
experimental and simulation curves follow the Arrhcnius 
behavior predicted by Eq. (0). In these plots the slope 
of the curve corresponds to the activation energy E c tf 
for the rate limiting step of the nucleation process. The 
results are E cS = 0.108 ± 0.005 and 0.112 ± 0.008 eV for 
the experimental and simulation curves, respectively. 

In Fig. ||we present the simulation results for the island 
density N vs. the coverage 9, for F = 4.8 x 10~ 4 ML-s -1 
and different temperatures within the experimental range 
213 - 263 K. The data extracted from the simulations is 
marked by *, while the connecting lines are only guides 
to the eye. The quick rise of N in the nucleation regime, 



T(K) 


Di (sites s ) 


Di/Lh. 


D s /D 1 


213 


3.6 ± 0.1 


1.5 ± 0.1 


0.021 ± 0.002 


223 


11.7 ± 0.4 


1.3 ± 0.1 


0.023 ± 0.002 


263 


520.0 ±15.0 


1.27 ±0.08 


0.040 ± 0.003 



TABLE II. The adatom diffusion coefficient and the rela- 
tive dimer and trimer diffusion coefficients, for the three stud- 
ied temperatures. Note that through the whole temperature 
range the relative diffusion coefficients almost do not change. 
The diffusion coefficients of islands of size 4 and higher are 
orders of magnitude lower. 

is followed by a plateau in the aggregation regime, and 
then a slow decrease that is the mark of the coales- 
cence regime. The coalescence at the lower temperatures 
T = 213 and 223 K, is system-size independent and is 
due to island separation being comparable to island size. 
At the higher temperatures, especially T = 263 K, the 
coalescence is mainly due to island separation being com- 
parable to system size, i.e., finite size effects. We be- 
lieve that this effect reproduces to a certain extent the 
effect of the steps limiting the width of the terraces in the 
experimental system.113 Our simulations indicate that at 
9 = 0.3 ML, where the experiments were done, the sys- 
tem is already in the coalescence regime. Therefore, the 
comparison of both experiment and simulations with the- 
ory [Eq. (|l|)] is misleading. 

The comparison between the experimental and simu- 
lation results shows that the set of energy barriers used 
here, provides a good quantitative description of the sub- 
monolayer growth in the low coverage regime. We now 
proceed to analyze the results in the context of available 
theoretical work. In Table |n| we present the calculated 
diffusion coefficients of adatoms as well as of small islands 
(dimers and trimers) at the experimentally relevant tem- 
peratures. The diffusion coefficient of a dimer D2, is 
found to be comparable to the adatom diffusion coeffi- 
cient D\. The diffusion coefficients of trimers and larger 
islands (not shown) are almost two orders of magnitude 
smaller. Thus, only adatoms and dimers are mobile on 
the surface on a time scale relevant to the growth process. 
In this case the meanialand separation in the aggregation 
regime is given by:E3E3E3 

e~( Dl D 2 /F 2 ) 1/W . (3) 

In Ref. ^ it was shown that the mobility of monomers 
and dimers is fully determined by only three energy bar- 
riers, the single adatom hopping Eq, the dimer lateral- 
bond breaking E2 and the re-establishing of a NN bond 
E± (the labeling follows Refs. 

£0 = ^(0) £ 2 =£(^) E^E^y 

To determine also the dimer stability one needs to specify 
the dimer bond-breaking energy Eg , and the next-nearest 
neighbor bond-breaking energy E\ : 
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Specifically, in the case of Cu, E = 0.485, E 1 = 0.563 
(In Ref. [16| the value of this barrier was mistyped.), 
E 2 = 0.463, E 4 = 0.183 and E 8 = 0.811 eV. Dimer break- 
ing can take place either from the NN configuration with 
energy barrier Eg, or from the next-nearest neighbor con- 
figuration with effective energy barrier E% — E4 + E\ = 
0.843 eV. The high values of these barriers guarantee the 
dimer stability. The adatom diffusion is a single step pro- 
cess, hence, Z?i = z/exp (~E q /KbT) [Eq. (@)]. The dimer 
diffusion is a double-step process involving a move with 
barrier E2 , followed by a move with E4 . The former move 
is the rate-limiting move, thus, up to a combinatorial fac- 
tor of order 1 we approximate D2 — vexp (—E2/ K-qT). 
Since E 2 — Eq we expect D2 and D% to be comparable, 
as is indeed the case (see Table [n]). Note that this en- 
ergy structure is predicted to be common to most of the 
fcc(001) metal surfaces.EEl 

Inserting D\ and D<x into Eq. (|3|) and comparing with 
Eq. (|) we get 

E cS = {E + E 2 )/10 = 0.0948eV, (4) 

that differs by more than 10% from the values obtained 
from the simulations (E e g — 0.112 ± 0.008 eV) and the 
experiment (E c g = 0.108 ±0.005 eV) presented above in 
Fig. ||. We recall that the evaluation of E e g was done at 
9 = 0.3 ML where the system is already in the coalescence 
regime, and therefore deviations from Eqs. (|l|),(||) and (^) 
are expected. Therefore, to examine the scaling relations 
of Eqs. (|J) and (||) we recalculate the island separation £ 
vs. 1/T in the same temperature range for a lower cover- 
age of 9 — 0.125 ML where the validity of those equations 
is better satisfied. We find that E cS = 0.096 ± 0.003 eV 
in good agreement with Eq. (^). These results are in 
excellent aeffieement with recent calculations by Boisvert 
and Lewis. c3 

We will now revisit the calculation of 7. According to 
the scaling theory [Eq. (||)] ^=2/5 for a system with mo- 
bile monomers and dimers.E3o The experimental results 
for 7 and the corresponding simulation results, obtained 
at 9 = 0.3 ML and presented in Table | are significantly 
different from that value. The discrepancy is removed 
when we recalculate 7 for T = 263 K and a much lower 
coverage, 9 = 0.125 ML, where Eq. (||) applies. The re- 
sults of this recalculation are presented in the rightmost 
column of Tabic |. Indeed we find 7 ~ 2/5, in good 
agreement with the scaling prediction. For T — 213 and 
223 K in the studied flux range the aggregation regime 
is wiped out, and the system is dominated by either nu- 
cleation or coalescence. Therefore, deviations from the 
scaling prediction are expected even for the lower cover- 
age 9 = 0.125 ML, as observed in Table |. 

The diffusion and coarsening of large Cu islands on 
Cu(001) has recently been studied experimentally.Efl The 
results were analyzed using MC simulation and a model 
that is qualitatively consistent with our model for the 



barriers relevant to large island diffusion J2a However, the 
models differ significantly in the lateral-bond breaking 
barrier, E 2l that is dominant in the monolayer growth 
through island nuclcation. Actually, the crucial factor is 
the difference E 2 - E = 0.463 - 0.485 = -0.022 eV in 
our model, and 0.52,— 0.399 = 0.121 eV in the model 
of Heinonen et al.E3 It follows that while D2/D1 = 
cxp[-(£; 2 - E )/K B T] is between 3.3 and 2.6 in the 
present model for temperatures between 213 and 263 K, 
it is between 0.001 and 0.005 in the model of Heinonen 
et al£3 Thus, according to their model dimers are static 
entities on the surface giving rise to a different scaling 
properties. On the other hand, the barriers of Shi et 
al.Jla E = 0.503 and E 2 = 0.494 eV, are in good agree- 
ment with ours. 

In summary, we have presented a comparison between 
kinetic MC simulations and experiments of Cu/Cu(001) 
growth in the submonolayer regime. The available ex- 
perimental data is for 9 = 0.3 ML and T = 213, 223 
and 263 K. At this coverage a good agreement between 
simulation and experimental results is found. However, 
these results deviate from mean-field predictions, a fact 
that we attribute to the system being away from the ag- 
gregation regime at 9 = 0.3 ML. Indeed, repeating the 
analysis of the simulation results for 9 = 0.125 ML we 
find 7 - 0.4 and E eS = 0.096 ~ (E + E 2 )/10 eV con- 
sistent with mean-field predictions when only adatoms 
and dimers are mobile on the surface. We conclude that 
for Cu/Cu(001), in the studied temperature range, only 
adatoms and dimers are mobile. Furthermore, mean-field 
predictions apply at 9 ~ 0.1 ML but not at much higher 
coverage. Finally, the comparison between simulation 
and experiment, augmented with a scaling procedure, en- 
ables us to depict the attempt rate v = 1.2 • 10 13 s _1 of 
the experimental system. 

We thank Hanoch Mehl for helpful discussions. Ex- 
perimental portions of this work were supported by Oak 
Ridge National Laboratory, which is managed by Lock- 
heed Martin Energy Research Corporation for the U.S. 
Department of Energy under Contract No. DE-AC05- 
960R22464. 
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FIG. 1. Comparison of experimental results (full symbols) and simulation results (empty symbols) for the island separation 
£ vs. the inverse flux 1/F, at three temperatures T = 213, 223 and 263 K. The coverage is = 0.3 ML. The solid lines represent 
fits to Eq. ([j]). The resultant values of the exponent 7 are presented in Table [l| 
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FIG. 2. The island separation I vs. 1/T, obtained from the experiment for F = 3.21 x 10 4 ML-s 1 and simulation for 
F — 4.8 x 1CP 4 ML-s -1 at coverage 9 = 0.3 ML. We also present the simulation results of the same runs but at an earlier 
stage when 9 = 0.125 ML. The solid lines represent fits to i ~ exp (-E ef t -/2K B T), with E eS = 0.108 ± 0.005, 0.112 ± 0.008 and 
0.096 ± 0.003 eV for the experiment, and the simulations at 9 = 0.3 and 0.125 ML, respectively. 
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FIG. 3. The island density N, vs. the coverage 9 at different temperatures with a flux F = 4.8 x 1CF 4 ML-s -1 . The solid 
lines are a guide to the eye. Since all the runs begin with an empty surface, the point (0, 0) is common to all the curves although 
it is not drawn in the figure. 



